How to formulate membrane potential in a 
spatially homogeneous myocyte model? 



A.J. Tanskanen^'^ , E.I. Tanskanen^ , J.L. Greenstein ^'^ and 

R.L. Winslow^ ^ 

^ The Center for Cardiovascular Bioinformatics and Modeling, The Johns 
Hopkins University School of Medicine and Whiting School of Engineering 
Baltimore MD 21218, USA 

^ The Whitaker Biomedical Engineering Institute, The Johns Hopkins 
University School of Medicine and Whiting School of Engineering, Baltimore, 

MD 21218, USA 

^Extraterrestrial Physics, NASA/Goddard Space Flight Center, Greenhelt, MD 

20770, USA 

Abstract 

Membrane potential in a mathematical model of a cardiac myocyte can 
be formulated in different ways. Assuming a spatially homogeneous my- 
ocyte that is strictly charge-conservative and electroneutral as a whole, 
two methods will be compared: (1) the differential formulation dV/dt = 
—I /Cm of membrane potential used traditionally; and (2) the capacitor 
formulation, where membrane potential is defined algebraically by the ca- 
pacitor equation V — Q/Cm- We examine the relationship between the 
formulations, assumptions under which each formulation is consistent, and 
show that the capacitor formulation provides a transparent, physically re- 
alistic formulation of membrane potential, whereas use of the differential 
formulation may introduce unintended and undesirable behavior, such as 
monotonic drift of concentrations. We prove that the drift of concentra- 
tions in the differential formulation arises as a compensation for failure 
to assign all currents in concentrations. As an example of these con- 
siderations, we present an electroneutral, explicitly charge-conservative 
formulation of Winslow et al. model (1999), and extend it to describe 
membrane potentials between intracellular compartments. 



1 Introduction 

The single most important variable in an electrophysiological whole-cell model 
is membrane potential, defined as the potential difference across the cell mem- 
brane. It drives both gating of ion channels and fluxes of ions through the 
membrane. The basis for electrophysiological single-cell modeling was first pro- 
vided by the Hodgkin-Huxley model of the neuronal action potential (1952). In 
their model, membrane potential was defined by an ordinary differential equa- 
tion in which the time-derivative of membrane potential equals the sum of all ion 
currents through the membrane divided by membrane capacitance. Their model 
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assumed constant ionic concentrations and dynamic concentrations were added 
in later models (DiFrancesco and Noble, 1985). While tracking of concentrations 
made the models more reaHstic, doing so introduced new problems, including 
drift of concentrations under repeated stimuH (see Fig. QJ, over-determined 
initial conditions and an infinite number of steady states (Guan et al., 1997; 
Hund et al., 2001; Kneller et al., 2002). 

Long-term drift of concentrations is present in the DiFrancesco-Noble (Guan et 
al., 1997) and Luo-Rudy models (Luo and Rudy, 1994). Using numerical sim- 
ulations, Hund et al. (2001) showed that the Luo-Rudy model has a steady 
state and no concentration drift, when it is assumed that the externally ap- 
plied stimulus current is carried by potassium (K^) ions and that this ion flux 
contributes to the rate of change in intracellular concentration. Similarly, 
Kneller et al. (2002) demonstrated that this flnding also holds true in their 
model. However, neither the reason for nor the origin of ion concentration drift 
has been addressed quantitatively. In this study, we explain the mechanism of 
concentration drift (Sec. K^.l|l and examine the number of steady states present 
in the above models (Sec. 

We also investigate alternative formulations of membrane potential which 
preserve charge-conservation and electroneutrality. In particular, we compare 
what we refer to as the differential and capacitor formulations of membrane 
potential (Sec. ISJ, and consider issues affected by the formulation of membrane 
potential (Sec. EJ. Two examples (Sections 01 and [^l show how membrane 
potential may be formulated rigorously in a computational model of the cardiac 
ventricular myocyte developed by Winslow et al. (1999; abbreviated WR.I). Our 
results show that the capacitor formulation provides a transparent, well-defined 
formulation of membrane potential in a spatially homogeneous myocyte model. 



2 The myocyte as a capacitor 
2.1 Formulation of membrane potential 

Following Hodgkin and Huxley (1952), a spatially homogeneous single-cell model 
is described as a parallel RC-circuit model. Given an initial value, membrane 
potential V is then defined by the differential equation 

- = -— /, (1) 

dt Cm 

where / is the sum of all outwardly directed membrane currents and Cm is total 
membrane capacitance. This formulation of membrane potential is referred to 
as the differential formulation. 

Membrane potential can be defined by assuming that a cell is a capacitor, 
consisting of an ideal conductor representing the interior of the cell, a dielectric 
representing the ability of the cell membrane to separate charge and a second 
ideal conductor representing the extracellular space surrounding the cell. In this 
capacitor formulation, membrane potential is defined as the ratio of charge Q 



2 



to total membrane capacitance C, 



V - Q/Cra- 



(2) 



Charge Q includes contribution from all charged particles in a cell. More specif- 
ically, membrane potential of a model containing a single intracellular compart- 
ment with N ion species {Sk}^^i is given by 



where zj- is the valence of species Sk, fmyo is volume of myoplasm and F is 
Faraday's constant. 

A biophysically detailed model often describes a myocyte using more than 
two compartments. The capacitor formulation ^ is extended to multiple com- 
partments by describing the cell as a network of capacitors (Fig. EK) ) where each 
dielectric corresponds to an interface between two compartments. The interior 
of a compartment is considered to be an ideal conductor insulated from other 
compartments by a membrane. The membrane potentials can be expressed 
as functions of charges of compartments and of capacitances of interfaces. In 
particular, a compartment that is completely enclosed within a larger compart- 
ment influences the membrane potential of the surrounding compartment only 
via charge, geometry is irrelevant (Griffiths, 1989). For example. Figures 
and |2j3 show graphical and capacitor representations of a three-compartment 
cell model. Membrane potentials are given by Vi — (Qmyo + Qmito)/Ci, and 
V2 = Qmito/C2, each depending only on net charge of the enclosed volume. 

2.2 Electroneutrality and homogeneous concentrations 

Membrane potential is generated by a small number of ions, the bulk ionic 
concentration is electroneutral (Hille, 2001). We require that (1) a model is 
electroneutral as a whole, that is, the net charge is zero; and that (2) the con- 
centrations are spatially homogeneous (i.e. well-mixed) in each compartment. 
Since such requirements are not necessarily satisfied in the presence of charged 
particles, we will show this to be the case for the model presented here. 

In an ideal conductor, induced charges balance electric field in such a way 
that potential is constant inside the conductor. Hence the interior of each com- 
partment is exactly electroneutral in the capacitor approximation and all net 
charge is located on the compartment boundary. This is consistent with nar- 
row physiological range of membrane potentials, which also requires that sum 
^Zk[Sk]i must be nearly zero. 

Both anions and cations must be present in each compartment, since the 
absence of anions (cations) would imply that all cations (anions) are located on 
the membrane. The concentration of ions that comprise the induced charge is 
extremely low compared to bulk concentrations and can typically be neglected. 
For example, in the WRJ model, range [—100,50] mV limits induced charge 




(3) 



k=l 
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to correspond to a monovalent ion concentration in range [—0.5, 0.25] /imol/L, 
which is neghgible compared to most intracellular ion concentrations. Elec- 
troneutrality of the interior allows the assumption that ion concentrations are 
homogeneous in each compartment. 

In a typical experimental setup, a myocyte is embedded in an electroneutral 
solution of essentially infinite volume relative to the cell volume. While this 
extracellular space (ES) is treated like any other compartment in a cell model, 
the assumption of infinite volume (implied by constant extracellular concentra- 
tions) brings certain extra complications that will be addressed in the following. 
If the myoplasm has charge Q, charge of the ES is —Q. However, if the charge 
density of ES, ^X^fcLi Zk[Sk\e (notation as in ijSl), is non-zero, then the ES has 
infinite charge which implies infinite membrane potential. Hence, the ES must 
have zero charge density, but typically non-zero charge. The paradox arises 
from infinite volume: if the charge of the ES is finite, the average charge density 
is zero because no finite fiow of ions can change the concentrations within an 
infinite volume eventhough charge is altered. As with all other compartments, 
all net charge is located on the boundary of the ES, and the remainder of the 
compartment is exactly electroneutral. In conclusion, the capacitor approxima- 
tion of the cell is consistent with the requirement of electroneutrality and the 
assumption that ions are homogeneously distributed in a compartment. 



2.3 Connection between the formulations 

Next, we will study the exact connection between the two formulations of mem- 
brane potential. In the differential formulation, membrane potential is defined as 
an independent variable and does not depend on ion concentrations, whereas in 
the capacitor formulation, membrane potential is a function of concentrations. 
Hence, in the differential formulation initial conditions are over-determined, 
since initial conditions are assigned independently for interdependent variables. 
This issue can be resolved by introducing implicitly-defined ion concentrations 
in the differential formulation, as will be shown in the following. 

For example, assume a one intracellular compartment model with M con- 
centrations {[S]m}m=i and N currents {Ik}k=i^ is defined in the differential 
formulation by 

^■=1 (4) 

1=1 

where v is the volume of the compartment, am.iJ = 1, ^m, indexes the Nm 
currents assigned to ion species Sm of valence z^. 

Assume that not all currents carry ion species modeled by any of the intra- 
cellular concentrations [S]„i- Such current are assigned to an implicitly-defined 
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monovalent concentration [S]i, 

-[S],^-Y,Ia,jFv, (5) 
1=1 

where Nj = N — J^m Combining equations 1^1 and © yields 

l^ = ^KEl(^4Sk) + |[Sh)/C™. (6) 

ni 

Integrating equation © gives membrane potential 

V = Fv{Y^ z^[SU + [S]i)/C,„. (7) 

m 

The initial value of [S]i is determined by the assumption regarding initial charge 
density. Equation 10 shows the exact connection of the differential formulation 
1^1 to the capacitor formulation J^J (in the two-compartment case, but the 
derivation can be generalized to any number of compartments). The major 
difference between the formulations is [S]i, which accounts for both the currents 
missing from equation Q and the charge missing from initial conditions in the 
differential formulation. The formulations are equivalent if [S]i is constant and 
all movement of charge is captured by the currents. Concentration [S]i and its 
time evolution are the source of a variety of problems, as will be shown in the 
following. 



3 Consequences of membrane potential formula- 
tion 

In this section, we will show how the physical constraints presented above in- 
fluence the formulation of a cell model. 

3.1 Concentrations and drift 

Here we examine how and under what conditions spurious concentration drift 
arises in the differential formulation. Assume a model with K"*" concentration 
[K]i and a constant anion concentration [B];, in which membrane potential is 
given by 

i^^'myo([K]i - [B]i)/C™. (8) 

It is clear that stimulus current /gtim must influence [K]i if it is to modify 
membrane potential. Charge-conservation and electroneutrality require that 
the stimulus current originates from the ES. 

In the differential formulation, stimulus current has not historically been 
assigned to a particular ion concetration. Formulating the above model in this 
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manner yields a model comparable to the one defined above defined by 



dV _ 1 

-TT — -J^K^K + -istimj, 
rf[K]i 1 , ^""^ 

~W " "f^^ — 

where /_r- is current through the cell membrane. Note that concentration [B]; 
does not appear in JHl. If we assign the stimulus current to implicitly-defined 
anion concentration [E]i (no information on valence of species E is contained in 
equation Q), 

and equation 10 implies that 

V^ = i^«myo([K]i-[E]i)/C„. (11) 

Since all currents are explicitly accounted for in the concentrations, equation 
lfTT|l shows that the model defined by equation is charge-conservative, when 
concentration [E]i is included in the formulation. A positive stimulus current 
decreases [E]i monotonically and eventually makes it negative, that is, [E]i has 
spurious drift. To keep membrane potential in the physiological range, the 
decrease of [E]i must be compensated for by a decrease of [K]\. Drift occurs 
in both the capacitor (equation lllll and differential formulations (equation |2l 
of the model, however, the cause of drift is transparent only in the capacitor 
formulation. 

Explanation of the drift presented above may be generalized to any number 
of concentrations. In a model with unassigned current Jc, membrane potential 
is given by equation The current Ic either increases or decreases [S]i mono- 
tonically. To keep membrane potential in the physiological range, charge density 
J2m^i ^"li^]™^ must compensate for the change in [S]i, which is manifested as 
a monotonic drift of at least one of the concentrations [S]ni. Even a minor drift 
in a single ion concentration will infiuence membrane potential if it goes un- 
compensated. Hence drift is inevitable if there exists a current that transports 
charge and hence affects membrane potential, but does not contribute to any 
ion concentration. This is consistent with previous numerical studies (Hund et 
al., 2001; Kneller et al., 2002). Previously presented descriptions of the drift 
(e.g., Hund et al. (2001)) have stated that lack of charge conservation is at the 
root of the problem but did not demonstrate how drift arises. As shown above, 
concentration drift is a compensation for the change in the implicitly-defined 
concentration [S]i. 

How does the drift depend on pacing rate? Since drift is relative to the 
charge transported by the stimulus current, higher pacing rate leads to more 
rapid decrease of [E]i, and consequently to faster drift in [K]i, consistent with 
numerical simulations of Hund et al. (2001). 



6 



3.2 Steady state and time scales 

An infinite number of steady states were observed in simulations of Guan et 
al. (1997), when initial concentrations were varied. If voltage is kept constant 
but concentrations are changed in the differential formulation, concentration 
[S]i (equation 10) is changed (Sec. |S3J- K [S]i remains constant during a 
simulation, different values of [S]i result in different steady states. Hence, an 
infinite possible number of initial values of [S]i yield an infinite number of steady 
states. On the other hand, [S]i is not constant if any of the currents implicitly 
assigned to [S]i is non-zero, in which case the model may not have a steady state. 
In the capacitor formulation, [S]i is explicitly included in the initial conditions, 
which resolves this issue. 

The time scale of exponential relaxation to steady state can be described by 
time constant r, measured by fitting A + Be"*/"^ to either the diastolic voltage, 
Na"*" or concentrations (Fig. 3). The time constant is mainly determined 
by the balance of Na"'' and K"'' concentrations. This time constant divides the 
model time scales into two regimes: slow (i > t) and fast (i < r). The fast times 
scale is on the order of a few seconds or less, whereas the slow time scale is on the 
order of tens or hundreds of seconds. In a study of fast time scale processes, Na"*" 
and K"*" concentrations can be clamped. However, if concentrations are clamped 
for a time period comparable to r, incorrect behavior will occur since implicit 
concentrations will change significantly during the simulation. Analogously, in 
voltage clamp a control current holds the membrane potential at the desired 
level. If this current is not assigned to any concentration, it is accounted for 
in implicit concentration and the behavior of the system will be incorrect in 
protocols with duration comparable to the time constant. 

3.3 Rapid equilibrium approximation 

The rapid equilibrium approximation (REA) often provides a powerful method 
to simplify a cell model (see, e.g., Hinch et al. (2004)). However, the REA 
requires instantaneous movement of charge in the model. Due to currents as- 
sociated with instantaneous movement of charge, the differential and capacitor 
formulations may yield different results, as will be shown in the following. 

Assume a model consists of a small-volume intracellular compartment and 
an infinite-volume extracellular compartment which are connected by a single 
ion channel that is described as a two-state Markov model. The intracellular 
compartment has anion concentration [A]i equal to extracellular concentra- 
tion [K]o. Also assume that the intracellular K"^ concentration [K]i is an instant 
function of the state of the channel (essentially a REA) : if the channel is closed 
[K]i = 0, otherwise [K]i = [K]o. Membrane potential V = F{[K\i - [K]o)/C™ 
depends on state of the ion channel. Assuming that a and (3 are opening and 
closing rates, respectively, of the ion channel, the probability Pc that the ion 
channel is closed evolves according to equation dPc/dt — (3 — (a + (3)Pc- 

Assume a large ensemble of N intracellular compartments described above to 
be located within a single cell. According to capacitor formulation, membrane 
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potential of the cell is then V{t) = —NFPc{t)[K]o/Cm- Since each individual 
channel allows non-zero flux only at the time of channel closing and opening, the 
total current seems to mostly be zero. However, the ensemble current through 
the cell membrane due to REA is / = F{(3 — (a + /3)Pc)[K\oN. This current is 
missed in the naive appHcation of the differential formulation, which thus yield 
incorrect membrane potential. 

In a recent reformulation of the WRJ model by Greenstein et al. (2005, sub- 
mitted to Biophys. J.), the REA was employed. In this case, current due to the 
REA is not easily solved from the equations, and use of the differential formu- 
lation introduces systematic, albeit small, error in calculation of the membrane 
potential (Fig. ^D-E). The small magnitude of this error is due to the small ra- 
tio of the diadic volume to the total cell volume. When a similar approximation 
is applied to a larger compartment, the difference can be signiflcant. 

4 Case study 1: the WRJ model 

The WRJ model describes intracellular Ca^^ dynamics of the cardiac ventricular 
myocyte. The model consists of four intracellular compartments: myoplasm, 
network sarcoplasmic reticulum (NSR), junctional sarcoplasmic reticulum (JSR) 
and the diadic space (SS). The myocyte is embedded in the extracellular space 
containing constant ion concentrations. The WRJ model has concentration drift 
and consequently does not have a steady state; it does not exhibit physiologically 
realistic dependence of APD on pacing rate, nor is it electroneutral. As an 
application of the considerations in the previous two sections, we reformulate 
the WRJ model and address these issues in the following. 

To address the above mentioned issues, we modified the model of (Mazhari et 
al., 2002) which was based on the WRJ model in the following manner: stimulus 
current /gtim is assigned to intracellular K"*" concentration [K]i; conductances of 
background calcium current /ca.b, background Na"*" current /Na,b and Na"*"- 
K"*" pump /NaK were adjusted to balance concentrations in the long term; and 
K^, Na"*" and stationary anion concentrations were added to compartments if 
not already present. We reformulated membrane potential using the capacitor 
formulation, in which membrane potential is given by 

V = [Q^ + Qnsr + QjSR + Qss]/C™, (12) 

where charges Qk are defined through 

Q^ = ^^«myo([K]i + [Na]i + 2[Ca]*°* - [S],), 
Qnsr = ^wnsr([K]nsr + [Na]NSR + 2[Ca]NSR - [S]nsr), 
QjSR = i^^^JSR([K]jsR + [Na]jsR + 2[Ca]*°4 - [S]jsr), 
Qss = J^«ss([K]ss + [Na]ss + 2[Ca]|°s* - [S]ss), 

where Vk is volume of compartment k, concentration [S]k is static anion con- 
centration in compartment k. The index tot refers to the sum of buffered and 
free Ca^+. 
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In a ventricular myocyte, action potential duration (APD) depends on pacing 
rate. Figure Q]^ shows steady state action potentials at pacing rates 0.25-2 Hz, 
with higher pacing rates producing shorter action potentials. The pacing rate 
dependence of APD arises mainly as a result of long-term changes in [Na]i and 
[K]i. Figure nj3 shows simulations started (1) with intracellular concentrations 
set to extracellular concentrations; (2) from a minimally perturbed steady state 
(static charges unchanged); and (3) with steady state [Na]i and [K]; replaced by 
[Na]i -I- 1 and [K]i — 1. In each case, limit cycles in {V, [Na]i) phase space are 
identical within numerical accuracy, showing that the steady state is unique, 
given parameters including pacing rate. 

Figure illustrates long-term changes in the model, sampled at 1 Hz. The 
simulation is started from a state with no concentration gradients. Figure 
shows diastolic membrane potential. Figure exhibits the increase of [Na]; 
(initially set to [Na]o) and the decrease of [K]i (initially set to [K]o) to their 
steady state values. Homeostasis is approached approximately exponentially 
with a time constant of 90 seconds. Figure shows the decrease of diastolic 
[Ca]i towards its steady state. Two oscillatory regimes of [Ca]i emerge: first 
at roughly 800 seconds with small, subthreshold oscillations; second with large 
amplitude oscillations at roughly 1,400 seconds. Figurel^t) shows the non-trivial 
time evolution of APD during the simulation. 

5 Case study 2: Compartmental membrane po- 
tentials 

Intracellular compartments are important for proper myocyte function. In par- 
ticular, the sarcoplasmic reticulum (SR) stores Ca^"*" ions. While the process of 
uploading of Ca^"*" into the SR is electrogenic, it does not seem to infiuence SR 
membrane potential that is observed to be small in amplitude (wsr; Bers, 2001). 
Pure diffusion of ions through SR membrane cannot explain ■^sr, since it would 
balance ion species separately without consideration to vsr- A small vsr re- 
quires that movement of counter-ions, likely and Cl~, balances the potential 
difference generated by Ca^^ movement (Pollock et al., 1998; Kargacin et al., 
2001). Indeed, SR membrane is known to have Cl~ channels gating according 
to Vsr (Townsend and Rosenberg, 1995), suggesting that vsr does affect Ca^"^ 
handling unless kept at almost zero voltage. Somewhat counter-intuitively, zero 
"^SR requires currents driven by vsr- 

5.1 Formulation 

To better understand the basis for vsr and the concentrations of ions in SR, we 
extend the model of Case study 1 to describe intracellular membrane potentials. 
In this case, the cell is described as a network of capacitors shown in Fig. |33. 
Membrane potential Vk of interface k (Fig. [33) can be expressed as a function 
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of capacitance Ck of and charge qk bound to interface k by 



(15) 



Vk = Qk/Ck. (14) 

Given net charges Qm = Fv„i J2i zi[Si]m, the charges qk are 

^ -P{b + C5S + Qjsr) - Qnsr/C? - s 
r + + C37 + Csr) 

92 = -qi - Qe: 

93 = 91C3/C1 - 92C3/C2, 

94 = 92 - 93 - QSS, 

95 = qaCs/C^ ~ q^C^/d, 

96 = Qjsr - 94 + 95, 

97 = -96 — Qnsr, 

where /3 = + 7 = c?7 + c^, & - Qss + Qe(l + gf ), r = [1 + 7(^^3 + 
C4)]/C4, s = {Qe{CiC2 + C2C4 + Ca) + C2Qss)l{CiC2). Charge Qe is given by 
electroneutraUty, Qe = -(Qm + Qjsr + Qnsr + Qss)- 

An ion channel senses voltage V across the membrane in which the channel 
is located. Then the flux Jq of ions of species S between compartments a and 
b is given by the Goldman equation (Hille, 2001) 

RT e^VF/RT _ I • (16) 

where D is the diffusion coefficient, z valence, T temperature and R universal 
gas constant. Ion channels (other than Ca^"^ channels) between intracellular 
compartments are assumed to always be open and to be selective to a single 
ion species. In particular, myoplasm-SS and JSR-NSR interfaces have large- 
conductance, permanently-open pores. Each compartment contains Cl~ and 
stationary anion concentrations. 

To include the effect of voltage V7 modulating Ca^"*" transport to NSR, we 
derived a model for the SERCA2a pump assuming Michaelis-Menten kinetics, 
that SERCA2a achieves equilibrium with Ca^"*" instantaneously, and that rates a 
and (3 out of the intermediate Michaelis-Menten state are related by ^-'^^'^1^'^ = 
a/P, where rjAG is free energy of ATP breakdown. The flux through SERCA2a 
is 



JSERCA - V^a. ^ ^ j^^j2/^2 + [Ca]^,3^/i^^,^ ' ^^^^ 



where Knsr = K,e^'^'^'^^^-^p-'^''°^-'^^'''^ , in which 77 = 0.74, cq is elementary 
charge, Ki = 300 nmol/L. 



5.2 Results 

The model has a realistic dependence of APD on pacing rate, as evidenced by 
Figure 0)4.. Intracellular membrane potentials are consistently small but non- 
zero (Fig. ^). Membrane potentials V3 and V7 between myoplasm and SS and 
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NSR are essentially zero. Membrane potentials V4, V5 and Vg between JSR and 
other compartments are non-zero, roughly 1 mV, as a result of Ca^"*" buffering 
by Calsequestrin in JSR, the absence of which would reduce these membrane 
potentials to essentially zero. External membrane potential of SS (V2) tracks 
closely external membrane potential of myoplasm, Vi. 

When diffusion rates between myoplasm and SR are reduced, intracellular 
membrane potentials are increased (Fig. EP)- The maximal NSR membrane 
potential that SERCA can overcome by employing energy of ATP hydrolysis 
(56 kJ/mol; Bers, 2001) is 237 mV. Without the movement of counter-ions, 
7.5 /imol/L would be the maximal increase of [Ca]NSR over [Ca]i, while the 
measured concentration difference is roughly 0.7 mmol/L (Bers, 2001). Since 
cells are sensitive to any disruption in Ca^^ handling, even a minor SR mem- 
brane potential has functional consequences. 

Case study 2 demonstrates that (1) Intracellular membrane potentials can 
be incorporated in a cell model, and that (2) they give physical constraints 
on intracellular concentrations and require delicate balance of charges; (3) the 
magnitude of SR membrane potential vsn can be explained by movement of 
counter-ions; and (4) buffering of Ca^^ affects JSR membrane potential. 

6 Discussion 

6.1 Which formulation is appropriate? 

The capacitor formulation expresses membrane potential as a function of charge 
and capacitance, whereas the differential formulation uncouples membrane po- 
tential from concentrations. The main differences between the two formulations 
are an integration constant (the differential formulation is time derivative of the 
capacitor equation), and "independence" of membrane potential from concen- 
trations in the differential formulation. These two issues imply the presence of a 
dynamic, implicitly-defined ion concentration(s) in the differential formulation 
(Sec. 13.111 . which makes interpreting simulation results difficult and prone to 
errors, in addition to the presence of spurious drift of concentrations and issues 
with steady state. The differential formulation is equivalent to a particular ca- 
pacitor formulation, which, however, may not be the one intended due to the 
presence of implicit concentrations. 

The capacitor formulation requires "almost" electroneutrality and carefully 
chosen initial conditions due to sensitivity of membrane potential to net charge. 
However, these are physical constraints, since even a small additional charge 
can have a drastic effect on membrane potential. The capacitor formulation 
provides a transparent formulation of membrane potential, and requires no 
implicitly-defined concentrations. The differential formulation is best viewed 
as a shorthand notation for the more complete and better-defined capacitor 
formulation. 
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6.2 Comparison with previous studies 

Guan et al. (1997) show that the DiR-ancesco-Nobel model has infinite num- 
ber of steady states, when initial concentrations are varied. To address this 
issue, they suggest that concentrations should be treated as parameters. How- 
ever, they neglected the implied concentration in their study, inclusion of which 
resolves the issue. Consistent with our explanation (Sec. Ki.lll . Kneller et al. 
(2002) observed that if the sum of concentration changes in initial conditions is 
zero, the steady state stays the same. 

The sinoatrial node model by Endresen et al. (2000) uses the capacitor equa- 
tion for membrane potential. However, the reasoning behind their definition is 
different from ours. The membrane potential of a simplified model described 
in Sec. 13.11 is given by equation i^Wmyo([K]i — [B]i)/Cm. Our interpretation of 
this equation is that concentration [B]; represents intracellular concentration of 
an anion concentration required for electroneutrality. In Endresen et al. (2000), 
^"''^myo[B]i/Cm represents charge of the extracellular space, which however can- 
not be computed directly from the concentrations in the extracellular space due 
to its infinite volume. 

Varghese and Sell (1997) derived a "latent conservation law" showing that 
membrane potential can be expressed as a function of concentrations. They 
showed that when each current is assigned to contribute to some concentration, 
the implicit concentration is constant. As they point out, the reason for emer- 
gence of the conservation law is best understood from the capacitor equation 
JSJ. However, they did not consider time-dependent changes in the implicit 
concentration, study of which is required to understand model phenomena such 
as concentration drift. 

Using numerical simulations Hund et al. (2001) showed that assignment of 
stimulus current to [K]i is sufficient to remove concentration drift in the Luo- 
Rudy model (Luo and Rudy, 1994). They interpret the drift as a consequence of 
"non-charge-conservative" formulation of the model. However, we proved that 
a "non-charge-conservative" formulation is actually charge-conservative (Sec. 
I2.3|l . when the implicit concentrations is taken into account, even in the pres- 
ence of spurious concentration drift (Sec. 13.111 . In particular, this means that 
the equations defining the differential formulation require that the system is 
always charge-conservative. On the other hand, a model can be "non-charge- 
conservation", which, however, does not necessarily indicate the presence of 
concentration drift. Hence, explanation of concentration drift requires study 
of charge densities implied by, but not explicitly included in, the differential 
formulation (Sec. I3.1|l . Hund et al. (2001) further state that the differential 
and capacitor formulations are always equivalent. However, that is true only if 
the currents present in the model capture all movement of charge, which is not 
necessarily the case in, e.g., rapid equihbrium approximation (Sec. 

In this article, we showed that the capacitor formulation provides a physically 
consistent, well-defined formulation of membrane potential, and it avoids the 
problems all too often found in the differential formulation. In conclusion, we 
see little reason to use the differential formulation as the definition of membrane 
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potential in a spatially homogeneous myocyte model. 



Acknowledgements 

The cell models are available at website of the Center for Cardiovascular Bioin- 
formatics and Modeling (http://www.ccbm.jhu.edu/). AT wishes to thank 
Dr. Reza Mazhari, Dr. Sonia Cortassa and Tabish Almas for helpful discus- 
sions. This study was supported by NIH (ROl HL60133, ROl HL61711, P50 
HL52307), the Falk Medical Trust, the Whitaker Foundation and IBM Cor- 
poration. The work of ET was funded by the National Research Council and 
Academy of Finland. 



References 

Bers, D.M., 2001. Excitation-Contraction Coupling and Cardiac Contractile 
Force, Second edition. Kluwer Academic Publishers. Dordrecht. 

DiFrancesco, D., Noble, D., 1985. A model of cardiac electrical activity incor- 
porating ionic pumps and concentration changes. Phil. Trans. R. Soc. 
Lond. B. Biol. Sci. 307:353-398. 

Endresen, L.P., Hall, K., H0ye, J.S., Myrheim, J., 2000. A theory for the 
membrane potential of living cells. Eur. Biophys. J. 29:90-103. 

Griffiths, DJ., 1989. Introduction to electrodynamics, second edition. Prentice- 
Hall, Englewood Cliffs NJ. 

Guan, S., Lu, Q., Huang, K. 1997. A discussion about the DiFi-ancesco-Noble 
model. J. Theor. Biol. 189(1) :27-32. 

Hille, B., 2001. Ion channels of excitable membranes, third edition. Sinauer, 
Sunderland MA. 

Hinch, R., Greenstein, J.L., Tanskanen, A.J., Xu, L., Winslow, R.L., 2004. A 
simplified local control model of calcium induced calcium release in cardiac 
ventricular myocytes. Biophys. J. 87(6):3723-36. 

Hodgkin, A.L., Huxley, A.F., 1952. A quantitative description of membrane 
current and its application to conduction and excitation in nerve. J. Phys- 
iol. (Lond.) 117, 500-544. 

Hund, T.J., Kucera, J. P., Otani, N.F., Rudy, Y., 2001. Ionic charge con- 
cervation and long-term steady state in Luo-Rudy dynamic cell model. 
Biophys. J. 81, 3324-3331. 

Kargacin, G.J., AH, Z., Zhang, S.J., Pollock, N.S., Kargacin, M.E., 2001. Io- 
dide and bromide inhibit Ca^"*" uptake by cardiac sarcoplasmic reticulum. 
Am. J. Physiol. Heart Circ. Physiol. 280(4) :H1624-34 



13 



Kneller, J., Ramirez, R.J., Charticr, D., Courtemanche, M., Nattcl, S., 2002. 
Time-dependent transients in an ionically based mathematical model of 
the canine atrial action potential. Am. J. Physiol. Heart Circ. Physiol. 
282: H1437-H1451. 

Luo, C.H., Rudy, Y., 1994. A dynamic model of the cardiac ventricular action 
potential. I. Simulations of ionic currents and concentration changes. Circ. 
Res. 74(6):1071-96. 

Mazhari, R., Greenstein, J.L., Winslow, R.L., Marban, E., Nuss, H.B., 2001. 
Molecular interactions between two long-QT syndrome gene products, 
HERG and KCNE2, rationalized by in vitro and in silico analysis. Circ. 
Res. 89(l):33-38. 

Pollock, N.S., Kargacin, M.E., Kargacin, G.J., 1998. Chloride channel blockers 
inhibit Ca^"*" uptake by smooth muscle sarcoplasmic reticulum. Biophys. 
J. 75:1759-1766. 

Townsend, C, Rosenberg, R.L., 1995. Characterization of a chloride chan- 
nel reconstituted from cardiac sarcoplasmic reticulum. J. Membr. Biol. 
147(2):121-36. 

Varghese, A., Sell, AR., 1997. A conservation principle and its effect on the 
formulation of Na-Ca exchanger current in the cardiac cells. J. Theor. 
Biol. 189, 33-40. 

Winslow, R.L., Rice, .J., .Tafri, S., Marban, E., O'Rourke B., 1999. Mechanisms 
of altered excitation-contraction coupling in canine tachycardia-induced 
heart failure, II: model studies. Circ. Res. 84(5):571-86. 



14 




Figure 1: (A) When stimulus current is assigned to K"*" concentration in Case 
Study 1 model, the model has a limit cycle ("steady state"; grey cycle) in 
([K]i, [Na]i) phase space. When the stimulus current is not assigned to any 
concentration, the model exhibits monotonic drifts away (black solid line) from 
the steady state; (B) Steady state action potentials at 0.25-2 Hz pacing rates in 
the Case Study 1 model. Membrane potential (ordinate; mV) is plotted against 
time (abscissa; ms); (C) Three simulations started from three different initial 
conditions (each with zero net charge) all approach the same steady state limit 
cycle in ([Na]i, V) phase space in the Case study 1 model. Sodium concentration 
(ordinate; mmol/L) is plotted against membrane potential (abscissa; mV); (D) 
Action potentials simulated using the differential (solid grey line) and capacitor 
formulations (dashed dark grey line) and (E) the potential difference between 
the formulations in Greenstein et al. (2005, submitted to Biophys. J.) model. 
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Figure 2: (^4) Spatial representation and (B) capacitor representation of a my- 
ocyte model with three compartments: extracellular space, myoplasm and mito- 
condria; (C) Capacitor representation of five-compartment single-cell myocyte 
model, consisting of myoplasm, subspace (SS), junctional sarcoplasmic reticu- 
lum ( JSR) , network sarcoplasmic reticulum (NSR) compartments and an extra- 
cellular space. Membrane potentials between the compartments are denoted by 
Vk. 
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Figure 3: A Case study 1 simulation run demonstrating how a myocyte makes 
its way from a state with no concentration gradients to physiological steady 
state: (A) Resting membrane potential (ordinate; mV) sampled once a second, 
plotted against time (abscissa; s); (B) Sodium (gray) and Potassium (black) 
concentrations (ordinate; mmol/L) against time (abscissa; ms); (C) Diastolic 
calcium concentration (ordinate; mmol/L; logarithmic scale) plotted against 
time (abscissa; ms); (D) Action potential duration (ordinate; ms; logarithmic 
scale) plotted against time (abscissa; ms). 



17 




Figure 4: Results from the Case study 2 model: (A) Shape of baseline action 
potential, paced at 1 Hz (solid black) and at 2 Hz (dashed grey). Membrane 
potential (ordinate; mV) is plotted against time (abscissa; ms); (B) Baseline 
membrane potentials (ordinate; mV) on myoplasm-JSR, myoplasm-NSR and 
myoplasm-SS interfaces plotted against time (abscissa; ms). (C) Increased 
membrane potentials (ordinate; mV) on myoplasm-JSR, myoplasm-NSR and 
myoplasm-SS interfaces plotted against time (abscissa; ms) in a simulation with 
reduced diffusion between intracellular compartments. 
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